Phase imaging from defocus information in a light field 
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Abstract 

Optical microscopy is without a doubt an essential component of life science research, but many objects 
of interest in biology are transparent. Chemical or immunological dyes, which can often be toxic, fluorescent 
transgenes, which require a protocol for transformation of exogenous DNA, and phase contrast, which 
inextricably combines phase information with amplitude, all serve to increase the contrast of transparent 
objects. Transparent objects of non-uniform refractive index and/or thickness are phase objects and phase 
can be extracted from defocus information above and below the object. The light microscopist will recognize 
this in that object contrast can vary as the stage is moved around an object's focal plane. This phase 
information can be quantitatively retrieved from a set of defocused images, but this has the disadvantage 
of decreasing temporal resolution, as the microscope stage has to be moved between image captures. By 
incorporating additional optics in the form of an array of pinholes or microlenses, the direction and position 
of all the rays can be recorded in a light field, a 4D function. By recording the light field, multiple planes of 
focus can be investigated, allowing for the capture of defocus information in real time. Here I discuss and 
investigate the use of the defocus information in a light field to obtain the transport of intensity equation, 
and use the fast Fourier transform method to solve this equation for phase. 

Introduction 



Biological samples of interest to the microscopist are 
often phase objects, that is, they are largely trans- 
parent and offer low contrast when viewed under 
hght illumination. Although transparent, they do 
perturb the phase of transmitted illumination via re- 
fractive index variation. This phase variation can be 
used to quaUtatively increase contrast as in Zernike 
phase contrast microscopy [5] through interferomet- 
ric means. In the past twenty years or so, phase has 
also been obtained quantitatively through a variety 
of means, one of the easiest being the transport of in- 
tensity equation's^ (TIE, equation[T]) which describes 
the change in intensity per change in axial distance 
from a given focal plane. The simplest way to ob- 
tain the transport of intensity equation is to record 
three or more images above, at, and below the fo- 
cal plane in which the object of interest resides [B]. 
Physically manipulating the microscope stage for de- 
focus decreases temporal resolution, so it is desirable 
to obtain defocus information in real time. Waller 
et al demonstrated that this can be accomplished by 



separating the color channels in a three color image, 
taking advantage of the chromatic aberration in mi- 
croscope objectives, specifically those that are either 
uncorrected for chromatic aberration or achromats, 
corrected for two colors [7j. Treating three color chan- 
nels as three separate images at slightly different fo- 
cus allows for the capture of quantitative phase in real 
time, allowing phase observation of dynamic events 
and video capture. The subject of study must be un- 
colored, making this method incompatible with say, 
fluorescence microscopy. Also, phase computation 
accuracy can be improved by incorporating many de- 
focused images [5], whereas recovery based on chro- 
matic aberration is limited to the color channels on 
a camera sensor. 

Light field microscopy can be achieved by inserting 
a microlens array with desired characteristics into 
the intermediate image plane of a microscope, 
and allows for the computation of multiple focal 
planes in a single capture [2 . This is strikingly 
amenable to real-time phase computation by solu- 
tion of the TIE, and has the advantage that many 
planes can be used to estimate the axial derivative 
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of irradiance. To demonstrate the capabihty to 
compute quantitative phase from defocus informa- 
tion contained in a Ught field, I generated a focal 
stack of images using a sample light field from the 
Stanford Light Field Microscope Project website 
fhttp:/ /graphics. stanford.edu/projects/lfmicroscope/) 
and built a MATLAB algorithm for solving for the 
associated TIE according to the fast Fourier 
transform methodas outlined in A, 



Method 

To demonstrate the feasibility of calculating phase 
information from defocus information contained in a 
light field, I utilized a sample light field available from 
the Stanford light field microscopy project jlj. Using 
the light field viewing software LFDisplay (available 
from the same source) I captured an image focused 
at zero as well as ten microns above and below uti- 
lizing screenshot capture. For experimental appli- 
cations outside of this concept note, data would be 
output directly from the light field manipulation soft- 
ware rather than as a copy of the monitor display. 

The intensity difference of the two defocused im- 
ages contains phase information retrievable via the 
transport of intensity equation[6]. The transport of 
intensity equationf4|P-222:e9.i2.2. 

^^«(^^_^ (1) 
oz 

Where ko is the wavenumber I{x,y) is the in- 



intensity on the right in Equation [T] is constant. This 

gives gP-222,e9.12.3. 



tensity. 



SIix,y) 
Sz 



is the difference in intensity for axial 



increments in z, and 4'{x, y) is the phase, ^^^g^^"^ is ob- 
tained empirically from a minimum of two defocused 
images [4]f ■222.-912.4, 
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Sz 2Az 
= 9{x,y) 



[I{x,y,Az)-I{x,y~Az)] (2) 



Where Az is the distance amount of defocus. To 
recover the phase from the difference image g{x, y) 
we assume, albeit counter-intuitively, that the image 
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V^c^{x,y) 



(3) 



To solve the Laplacian (V^) in Equation[31 we take 
the fast Fourier transform (FFT) of g{x,y) to yield 
G{kx,ky). In my MATLAB script, this is obtained 
with the 2-dimensional FFT command, fft2(). Phase 
in the frequency domain, ^{kx,ky) is obtained from 

the following [4jP-223,eg.l2.5. 



ko G{kx,ky) 



(4) 



Iq is the intensity distribution at the focal plane, 
assumed to be uniform. This value is obtained by 
taking the average of the in-focus image intensity. To 
convert the phase ^{x,y) in Equation |4] to the spatial 
domain, we take the 2-dimensional inverse FFT. 

To compute phase from the defocus stack obtained 
with LFDisplay I wrote a MATLAB function, which 
uses the fast Fourier transform to recover phase in- 
formation from the intensity derivative from a three 
image focus stack. The inputs are images above, be- 
low and at the focal plane of interest and the defocus 
distance between images, and the function returns a 
phase map. This function is included as an appendix 
to this note. 



Example 

A sample light field was used to test the feasibility 
of recovering phase. The subject was an onion skin 
taken with a 20X magnification, 0.5 numerical aper- 
ture objective, imaged onto a 24X36 mm array of 
square 125 micron pitch, f/20 microlenses, obtained 
from [1]. Three images separated by 10 microns of 
focus were used to calculate phase (Figure [1]). The 
difference image representing ^'^^g^^^ (amplified by a 
factor of 10 for viewing) and the recovered phase map 
are the contents of Figure [2j 
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Figure 1: Three 2-dimensional images separated by 10 microns of focus were 
retrieved from the Hght field. These will be used to compute the phase 

using the transport of intensity equation (TIE). Original LF image source: 
|http: / / graphics.stanford.edu/projects /Ifmicroscope / 2006.html 
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Conclusions 
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Figure 2: The defocus difference image (top, 
signal amplified lOX) is used to determine the 
function g(x,y) (Equation [3]), which is used to 
compute the phase map via Equation [4] (bot- 
tom). The fold in the onion skin is readily dis- 
cernible in the phase map, as well as certain 
details of the nodules in the vascular network 
of the plant. 



Any application of quantitative phase imaging could 
benefit by the use of light field imaging to push the 
technique into real time, for a relatively inexpensive 
cost of one or more interchangeable microlens arrays 
tailored to the application. The breadth of data col- 
lected in a light field could eliminate the requirements 
for an expert operator, for example in a rural clinic, 
because the phase computations and diagnoses can be 
made at a point spatially and/or temporally removed 
from the sample and microscope. Applications in 
manufacturing could include precise quality control 
of surfaces or structures differing only in thickness or 
refractive index, for example in laboratory etching of 
microfluidic devices in glass or silicon. 

There is a tradeoff in light field capture versus tra- 
ditional 2D microscopy, discussed in [2\ and [3]: the 
product of angular and lateral resolution cannot ex- 
ceed the diffraction limit of the optical apparatus and 
the wavelengths measured. In simple terms, the pixel 
area of each microlens image multipled by the total 
number of microlenses can not exceed the diffrac- 
tion limit of the microscope. This would limit the 
observable feature size to a few times larger than 
the resolution limit of a given optical arrangement, 
but may be attenuated with further computational 
tricks and should be a factor in microscope design, 
in particular in the choice of microlens array char- 
acteristics and objective. These considerations are 
discussed in a tutorial in [T]. I also note that Levoy 
et al have experimented with light field phase recov- 
ery of a homogeneous medium in [3 . However, their 
method employed additional optical components in a 
microlens array and led projector screen in the illu- 
mination path. 
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Appendix 



function [ phase ] = phaseNewWorld(Ia,IO, Ib,z) 

%generate a phase map based on the defocused images la (above), lb (below, 
%focused image 10, sind the defocused difference (in microns) z 

7, Detailed explanation goes here 

%For starters we will split of the green channel and deal with it 
"/.exclusively 

if size (size (la) ) == 3 °/o check for color/B&W image 

lag = Ia(l:size(Ia,l) ,l:size(Ia,2) ,2) ; 

Ibg = Ib(l:size(Ia,l) ,l:size(Ia,2) ,2) ; 

lOg = IO(l:size(Ia,l) ,l:size(Ia,2) ,2) ; 

else 

lag = la; 
Ibg = lb; 
lOg = 10; 
end 

y.Iavg is the uniform intensity distribution 
lavg = mean (mean (lOg) ) ; 

y.wavenumber kO (inverse microns, for green light) 
kOg = (2*pi())/(500 * 10--3); 

y.gofxy is the data that describe the change in field intensity per axial 
7odef ocus (z) . Ironically an assumption of this solution of the TIE is that 
"/othe field intensity is assumed to be constaint over z (see p. 223 Popescu 2011, 
y.QPI of cells and tissues) 
gof xy = (lag - Ibg) / (2*z) ; 

Xgofxy contains the phase information [gofxy = -lO/kO $\nabla"2$ \phi(x,y)] 
7oTo solve the Laplacian we convert gofxy to the frequency domain with the 
7.2D fft 

GOFkxky = fft2 (gofxy); 

7«debugging: show the frequency domain and difference image 

7,figure(l); subplot (2 , 1 , 1) ; title ('difference la-Ib'); imshow(Iag-Ibg) ; 

7.subplot(2,l,2) ; title ('FFT magnitude of dl/dz $\delta$'); imshow(abs (GOFkxky) ) ; 

y.Solve for the frequency domain phase: 

7oDeclare memory to hold the fft of phase 
PHI= zeros (size (lag, 1) ,size(Iag,2)) ; 
kx = 1; 

while kx < size (GOFkxky, 1) 
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ky = 1; 

while ky < size (GOFkxky , 2) 

PHKkx.ky) = (kOg/Iavg) * GOFkxky (kx , ky) / (kx"2+ky"2) ; 
ky = ky+1; 
end 

kx = kx+1; 

end 

phase = ifft2(PHI); 

figure(2); subplot (1 ,3 , 1) ; title ( 'defocused above ') ; imshow(Iag); 
subplot (1,3, 2) ; title( 'focused image ') ; imshow(IOg) ; 
subplot (1,3, 3) ; title ('defocused below') ; iinshow(Ibg) ; 
figure (3) ; 

subplot (1, 2, 1) ; title ('difference image amplified lOX, dl/dz'); imshow(10* (lag-Ibg) ) ; 
subplot (1,2, 2) ; title ('recovered phase, surf ace plot ') ; mesh(abs (phase) ) ; 
figure(4); title ( 'recovered phase, surf ace plot ') ; mesh(abs (phase) ) ; colormap(jet) ; 



end 
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